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Condensed matter physicists have long sought a realistic two-dimensional (2D) magnetic sys- 
tem whose ground state is a spin liquid — a zero temperature state in which quantum fluctuations 
have melted away any form of magnetic order. The nearest-neighbor S = | Heisenberg model on 
the kagome lattice has seemed an ideal candidate, but in recent years some approximate numerical 
approaches to it have yielded instead a valence bond crystal. We have used the density matrix renor- 
malization group to perform very accurate simulations on numerous cylinders with circumferences 
up to 12 lattice spacings, finding instead of the valence bond crystal a singlet-gapped spin liquid 
with substantially lower energy that appears to have Z2 topological order. Our results, through a 
combination of very low energy, short correlation lengths and corresponding small finite size effects, 
a new rigorous energy bound, and consistent behavior on many cylinders, provide strong evidence 
that the 2D ground state of this model is a gapped spin liquid. 

PACS numbers: 



I. INTRODUCTION 

A spin liquid (SL) is a magnetic system which has 
"melted" because of quantum fluctuations, even at zero 
temperature: 1 More precisely, it is a spin ground state 
with no broken symmetries, namely no preferred spin ori- 
entations, and no frozen valence bond patterns. Spin liq- 
uids, besides being exotic and novel forms of matter, have 
been suggested as possible routes to high temperature 
superconductivity,— and their unusual topological prop- 
erties may be useful in building quantum computers. 3 
A "Holy Grail" of condensed matter physics has been 
to find a realistic two-dimensional (2D) spin liquid. 
Spin-liquid-likc behavior is common in one dimension — 
quantum fluctuations are stronger in lower dimensions — 
and some systems with spatial anisotropy, in effect inter- 
polating between one and two dimensions, are thought 
to be spin liquids.— The strong charge fluctuations asso- 
ciated with being near the Mott transition, past which 
there are no localized spins, apparently can produce a 
spin liquid even without geometric frustration. 5 Strong 
multispin interactions, such as ring exchanges, are asso- 
ciated with this regime, and these terms have been in- 
cluded in 2D Hamiltonians to force spin liquid ground 
states^ Of great interest would be a spin liquid with a 
Hamiltonian that is isotropic in both space and spin di- 
rections, with no multispin interactions and preferably 
only nearest-neighbor exchange interactions on a simple 
lattice, where the spin liquid behavior is driven entirely 
by frustration. Here we present strong evidence for just 
such a spin liquid. 

The Heisenberg S = i kagome lattice model (HKLM) 
has long been thought to be the ideal candidate. On 
the experimental side, the material herbertsmithite, 
ZnCu3(OH)eCl2, is a spin-1/2 kagome antiferromagnet 
showing spin liquid behavior^— It is thought to be well 
modeled by a HKLM Hamiltonian with disorder and an 



additional Dzyaloshinski-Moriya interaction. The role of 
these perturbations in driving the spin liquid behavior 
is unclear»2r— On the theoretical side, a key problem is 
that there are no exact or nearly exact analytical or com- 
putational methods to solve infinite 2D quantum lattice 
systems. For one dimensional systems, the density ma- 
trix renormalization group 11 ' 12 (DMRG), the computa- 
tional method used in this work, serves in this capacity. 
Another "Holy Grail" has been to find an accurate and 
widely applicable computational method for 2D many- 
body quantum systems. 

The HKLM has been studied with approximate ap- 
proaches by leading theorists for decades^ with pro- 
posals for spin liquids and for valence bond crystals 
(VBCs)J^ In the last few years some new numerical evi- 
dence has suggested that the HKLM ground state is not 
a spin liquid. Instead, calculations have indicated that 
it might be a VBC, with a large 36-site unit cell, called 
the honeycomb valence bond crystal (HVBC). First pro- 
posed by Marston and Zengi^, and explored in more de- 
tail by Nikolic and Senthiii^, the HVBC was studied with 
a perturbative series expansion by Singh and Huse^, who 
found the series for the ground state energy to be rapidly 
converging and to give a low energy. Subsequently, Even- 
bly and Vidal utilized the multiscale entanglement renor- 
malization ansatz (MERA) on the HKLM to obtain a 
similar estimate of the energy of the HVBC 17 . MERA 
is related to DMRG but can be implemented directly in 
2D. This MERA study also produced the lowest numeri- 
cally exact upper bound to the energy per site, which was 
close to the series expansion energy. A recent DMRG 
study-i£, in contrast, found a spin liquid ground state, 
but on the largest lattices the energy obtained was sub- 
stantially above that of the HVBC, suggesting that the 
method had not found the true ground state. 

Here we perform a much more extensive DMRG study 
of the HKLM with important differences in techniques 



from the previous study. Most importantly, we avoid 
fully periodic (toroidal) boundary conditions, which are 
known to greatly magnify the truncation errors in DMRG 
compared to samples with open ends 12 . We present com- 
pelling evidence that the ground state of the HKLM is 
a spin liquid. This spin liquid has an energy which is 
substantially below that of the HVBC. We also provide 
a new rigorous upper bound on the ground state energy 
of the infinite 2D system. The spin liquid ground state is 
gapped to both singlet and triplet bulk excitations, with 
short correlation lengths. The short correlation lengths 
make the finite circumference (up to 12 bond lengths) 
cylinders that we can study converge rapidly towards two 
dimensions both in energy and properties. The spin liq- 
uid resembles a Z- 2 topologically-ordered short range res- 
onating valence bond (RVB) stated - — but has signifi- 
cant internal structure with a key role apparently played 
by eight-site resonant loops. 
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II. DMRG METHODS 

The density matrix renormalization groupii is ideal 
for one dimensional systems, but also has long been used 
for finite-width 2D strips or cylinders 23 ' 24 . To maintain 
a constant accuracy as the width of the system is in- 
creased, the number of states kept per block m must in- 
crease exponentially, increasing the computer resources 
needed, and thus limiting the maximum width feasible. 
High accuracy is required to distinguish between phases 
in the KLHM, so we have limited the systems studied 
to a maximum circumference of 12 lattice spacings, for 
which we can obtain a relative error of less than 0.1% in 
extrapolated energy. 

We use a cylindrical geometry for most of our calcu- 
lations, which we label by the cylinder's orientation and 
circumference. We put the circumference of the cylinder 
either along or close in orientation to the y-axis. For ex- 
ample, in our notation YC8 indicates a cylinder (C) in 
which one of the three bond orientations is along the y- 
axis (Y), and the circumference is 8 lattice spacings. For 
cylinders oriented with some bonds along the ^-direction 
(X), the circumference is measured in units of v3/2 times 
the lattice spacing, so that XC8 has a circumference of 
c = 4-\/3 lattice spacings. For cylinders which connect 
periodically around the cylinder's axis with a shift, the 
shift is added to the label: e.g. YC9-2 indicates a YC9 
cylinder which is connected with a horizontal shift at 
the "seam" of two units; YC9-2 has a circumference of 
c = 2v / 21. In these cases with shifts the kagome lattice 
is wrapped on the cylinder in a spiral. 

There are two key issues that complicate the calcula- 
tions: 1) the geometry of the cluster — the circumference, 
the periodic shift, the boundary conditions at the open 
ends — can affect the ground and excited states signif- 
icantly; and 2) on the wider systems, DMRG can get 
stuck in metastable states. Both of these problems can 
be surmounted, but it is necessary to study a wide variety 



FIG. 1: Valence bond pattern for the left two-thirds of a 208 
site XC8 cylinder. The width of the lines is proportional to 
KiSi-Sj)! (top two panels and left key) or \(Si-Sj)—e a \ (bottom 
panel, right key), with dashed lines indicating the quantity is 
positive. Here a signifies bond direction and the approximate 
bulk bond correlations are e a = —0.223 for horizontal bonds 
and —0.217 for diagonal. The top panel shows sweep 6 with 
m — 200, the middle panel sweep 14 with m = 600, and 
the bottom panel the final sweep 34 with m — 8000. Near 
m — 2400 (not shown), the energy drops fairly abruptly by 
about 0.1%, and the HVBC is replaced by a spin liquid state. 
The line widths have been constrained to a maximum near 
the left edge, bottom panel. 

of systems and initial states to find out which geometries 
do not frustrate the lowest energy states and how to pre- 
pare initial states which permit the algorithm to find the 
ground state. The results presented here are a very small 
fraction of the total number of systems studied, and have 
been selected to illustrate the key behavior most clearly. 



III. THE VALENCE BOND CRYSTAL VERSUS 
THE SPIN LIQUID 

In this section we present some of the key evidence 
that the ground state of the KLHM is a spin liquid, and 
specifically not the HVBC. To rule out the HVBC on 
the cylinders we have studied, despite the possibility of 
metastable states, a valuable technique is to favor the 
HVBC, both in initial state and boundary conditions. In 
spite of thus favoring it, the HVBC is unstable in the 
resulting simulations, which is strong evidence against 
an HVBC ground state. A typical example is shown in 
Fig. 1. In this simulation the cylinder circumference and 
wrapping vector accommodate the HVBC state, the left 
and right edges of the cluster were trimmed to accom- 
modate and pin the HVBC state, the initial state was 
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prepared in the HVBC state, and the ordering of the 
sites of the cylinder used by the DMRG followed an ir- 
regular path which always makes any two sites sharing a 
valence bond in the HVBC adjacent. (This biased spe- 
cial ordering allows a perfect non-resonating HVBC to be 
represented keeping only m = 2 states per block.) These 
biases towards the HVBC state (particularly the order- 
ing) make it metastable out to about m = 2400, after 
which it transitions to a spin liquid state. The irregular 
edges (along with a variety of other considerations) allow 
us to completely rule out the possibility that the uniform 
state is a superposition of shifted HVBCs. The final state 
is a remarkably uniform spin liquid in the center of the 
system, with only slight perturbations from the open, ir- 
regular edges and incomplete convergence, and a slight 
anisotropy from the finite circumference. This system 
clearly has a very small length scale for the decay of the 
perturbations due to the ends. If we use a more conven- 
tional ordering of sites, the HVBC is not even metastable; 
it disappears within the first few sweeps with m < 200. 
Energies are generally lower for a fixed m for the stan- 
dard ordering than for the special ordering. 

We do not find the HVBC to be the ground state on 
any of the cylinders we have studied. However, on the 
largest circumference cylinders, for which our accuracy is 
much reduced, the HVBC can be metastable even with 
the more unbiased standard site ordering. (A small bias 
is still present favoring the less entangled HVBC state.) 
For these systems we can compare energies between a pre- 
pared HVBC state and other states. Even in cases where 
the final energies are not accurate enough to clearly pick 
one state, we may be able to judge that the HVBC is 
not the ground state. In Fig. 2 we show such a case 
for a 400 site YC12 cylinder. Two runs were performed, 
one pinned to start as a HVBC, the other started with 
an essentially random initial state. Just after the pin- 
ning is released at m = 600, we see the HVBC has lower 
energy. However, near m — 2000 the initially random 
state is finding the spin liquid state and becomes lower 
in energy. Subsequently, the difference between the two 
energies widens as m is increased. The final state (not 
shown) appears to be a spin liquid with several localized 
defects which the simulation had not yet eliminated. 

The series expansions of Singh and Husei£ treat the 
HVBC strong-bond interactions as having a fixed J = 1, 
while the weak bonds have their interactions modified by 
a factor A; the expansion is about A = 0, and the uniform 
kagome lattice has A = 1. We have studied the ground 
state along this path of modified Hamiltonians H (A) for 
a 194-site XC8 cylinder. Along this path the system has 
an apparent first-order phase transition near A c = 0.984. 
This transition is most clearly seen using a "hysteresis 
plot". In Fig. 3 we show results from three different 
runs. Each run has a fixed m, and A is changed between 
each DMRG sweep, with A first increasing, then decreas- 
ing (note the arrows). The simulation cannot adapt the 
wavefunction from one phase to another in one sweep, so 
the system shows hysteresis, staying in the old higher en- 
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FIG. 2: Total energy E as the number of states kept m is 
increased for a 400 site YC12 cylinder, with two different ini- 
tial states — one essentially random, the other with temporary 
pinning fields to force an HVBC. 
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FIG. 3: Hysteresis plot for a 194 site XC8 cylinder show- 
ing the apparent first-order phase transition between valence 
bond crystal and spin liquid. For each of the three runs, 
the Hamiltonian parameter A (see text) is changed for every 
DMRG sweep, thus tracing out a path through parameter 
space that resembles the time evolution of a system under 
changing A. A linear approximation to the main trend of the 
energy E with A has been subtracted out. 



ergy state for several sweeps before drifting down to the 
new lower energy state. The crossing point of the curves 
in the two different directions converges very rapidly with 
increasing m to the transition point A c , even though the 
total energies have lower accuracy. This first-order tran- 
sition shows why the series expansion converged so welljiS 
but to an energy that was not the ground state energy 
at A = 1. 

For any of the cylinder geometries we can obtain an es- 
timate of the energy per site for the infinitely long cylin- 
der by subtracting energies of cylinders of different length 
to eliminate end effects. The best results are obtained by 
doing a sequence of m values for each cluster and extrap- 
olating in the truncation error. For A = 0.98 on XC8, 
which is in the HVBC phase, we obtain an energy per 
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site of —0.43237(4). This agrees fairly well with the se- 
ries expansion energy for this cylinder and A, —0.431(1). 
This supports the idea that the series expansion gives 
a reasonable estimate of the energy of the HVBC phase 
at A = 1 in two dimensions: —0.433(1) as does the 
MERA HVBC energy -0.4322^ which is a rigorous 
upper bound. MERA produces a rigorous upper bound 
because it generates a wavefunction for the infinite 2D 
system and evaluates its energy exactly (up to floating 
point round-off errors) 



IV. GROUND STATE ENERGIES 

It is possible to generate rigorous upper bounds on the 
ground state energy of the infinite 2D system from our 
results for finite open systems. Consider an open cluster 
C which can be "tiled" to fill all of 2D, with no sites left 
out, and having an even number of sites Nc- We take 
as a 2D variational ansatz a product wavefunction, the 
product being over all the tiles, where we use our best 
variational wavefunction for C (call it |C), with energy 
Ec) as the ansatz for each tile. The energy of any of 
the missing bonds connecting different tiles is zero, since 
(C\Si\C) = for any spin i. Therefore the energy per 
site of this simple product wavefunction is Ec /Nc ■ 

This approach is crude and converges slowly with the 
cluster size, with an error proportional to one over the 
width. Nevertheless, the SL energy is sufficiently low 
that we have been able to obtain a new rigorous upper 
bound on the 2D energy: E ( Q 2D) < -0.4332. This was 
obtained with a width- 12 open strip (which looks like 
XC12 unrolled) with Nc — 576, keeping m = 5000 states. 
The interior of this cluster had the uniform valence bond 
pattern expected for a spin liquid. 



TABLE I: Ground state energies and gaps for infinitely long 
cylinders of various circumferences, c. The third column 
indicates whether the diamond pattern fits perfectly on the 
cylinder. 
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FIG. 4: Comparison of energies per site for various lattices 
and methods. For cylinders, the horizontal axis in this plot 
is the inverse circumference in units of inverse lattice spac- 
ings. For tori ) 18 i 25 ~ ^£ the smallest circumference was used. 
In one case we show Lanczos energies for two different sized 
(36 and 42 sites) tori that have the same circumference . 26 ' 27 
The MERAii and our DMRG upper bound results apply di- 
rectly to an infinite two dimensional system, as does the series 
HVBC result that is plotted on the axis. The torus DMRG 
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energie: 

energies for those tori 



— are also upper bounds on the true ground state 



Our DMRG results are presented in Table I. The 
ground state energies are also plotted and compared to 
other calculations in Fig. 4. The DMRG energies are 
consistent with the Lanczos results^— and well below 
the energies of MERAil and the series expansions for the 
HVBCJ£ We note that the previous DMRG result is 
close to the true ground stated for the circumference 6 
torus. The entanglement across a cut that separates a 
circumference 6 torus into two parts should be roughly 
the same as across a cut that separates a circumference 12 
cylinder. We find that circumference 12 is presently our 
limit for obtaining good ground state energy estimates 
on cylinders. Thus it is perhaps not surprising that the 
DMRG results for torii£ give substantial overestimates of 
the ground state energies for circumferences larger than 
6. But these estimates may alternatively be viewed as 
variational upper bounds obtained with DMRG. 

The XC8 cylinder (1/c ~ 0.14) allows a direct com- 
parison of the energies between the HVBC series and our 
DMRG: the DMRG energy is lower by 0.004(1), and the 
series result for XC8 is near the 2D result. The corre- 
sponding torus shows much larger finite size effects in the 
HVBC series^ but the true finite size effects between 
the tori and cylinders arc quite small, as seen by the 
nearly identical results from Lanczos on tori and DMRG 
on cylinders when we use the largest available torus at 
each circumferencei^r— This is consistent with the small 
correlation length apparent in Fig. 1. We conclude that 
our widest cylinders would have minimal finite size ef- 
fects even if the system were in the HVBC phase; in the 



■5 



spin liquid phase, the finite size effects for the energy 
are much smaller, and our ground state energies for the 
largest cylinders apply to 2D with minimal correction. 

We now propose some nonrigorous bounds on the 
ground state energy of the infinite 2D system, by distin- 
guishing between frustrated and unfrustrated cylinders. 
We call a cylinder unfrustrated if its ground state en- 
ergy per site is lower than that of the 2D system. Here 
the periodic boundary conditions do not raise the ground 
state energy. Clearly, most of the very narrow cylin- 
ders have low energy (except possibly YC5-2), presum- 
ably due to the contribution to the energy from short 
RVB resonance paths that wind around the cylinder. We 
also believe that the "even" cylinders that are compat- 
ible with the 12-site unit cell diamond pattern that we 
discuss below are unfrustrated. Thus we propose that 
the energy of the XC8 cylinder provides a nonrigorous 
lower bound of —0.4382 per site on the 2D ground state 
energy. This cylinder has the same circumference as the 
36-site torus that is the "standard" exact diagonalization 
cluster (torus). This "XC8" torus has ground state en- 
ergy -0.43838 per site^ which is 0.00014(2) lower than 
the XC8 cylinder. The XC8 cylinder has only 1/3 as 
many extra winding resonant paths as does the torus, so 
should have a correspondingly smaller finite-size effect in 
its ground state energy. 

We believe that we also have at least one case (and 
probably many) of a cylinder that is frustrated, with a 
ground state energy above that for 2D. The argument 
for this also starts with the XC8 cylinder. Then we cut 
out a section of this cylinder with 42 sites and connect it 
to make a torus. This torus has various different new 
circumferences (topologically nontrivial ways to travel 
around it), which include YC7-2 and XC10-1. We sent 
this new (42-site!) torus out for exact diagonalization to 
Andreas Lauchli who promptly found its ground state en- 
ergy, which is approximately 0.0001 per site higher than 
that of the XC8 cylinder^ From this we conclude that 
one or both of the new circumferences YC7-2 and XC10- 
1 must be frustrating to the spin liquid ground state and 
thus that cylinder should have a higher energy than the 
2D system. This provides a nonrigorous upper bound on 
the 2D ground state energy of —0.4376. Thus our current 
estimate of the 2D energy per site is 

E^ 2D) = -0.4379(3) . 

Note that all the cylinders wider than YC8 have the same 
energy within the uncertainties, and thus rather small 
finite-size effects on their ground state energies. 

V. THE STRUCTURE OF THE SPIN LIQUID 
GROUND STATE 

The spin liquid ground state that we find has only 
short-range correlations and its ground state energy con- 
verges rapidly for finite circumference cylinders. We also 
find, as we discuss below, that there is a nonzero energy 




FIG. 5: Response of the spin liquid to small perturbations on 
cylinder YC8. In each case certain bonds have been strength- 
ened: upper left, an 8 site diamond by 1%; upper right, a 6 
site hexagon by 1%; and bottom, the wide strong and weak 
vertical bonds (28 of each) by ±0.5%. Line widths indicate 
subtracted bond correlations, as in Fig. 1, bottom panel. Sur- 
rounding dimers arise in response to the diamond pinning in 
the upper left panel; the diamonds arise in response to dimer 
pinning in the lower panel. 



gap for any excitations, including spin-singlet excitations. 
The structure of this ground state is apparently some sort 
of short-range resonating-valence-bond (RVB) stated 
The smallest resonant loops of singlet dimers in a nearest- 
neighbor RVB state on a kagome lattice each surround 
only one of the hexagons of the lattice. The possibilities 
are tabulated in, e.g., Table I of Nikolic and SenthilJ^ If 
the RVB state had all dimer covers equally weighted, all 
32 of these elementary loops would be equally present. 
What we find is that the ground state instead appears to 
substantially overweight certain 8-site loops of a diamond 
(or rhombus) shape. We first noticed this in seeing what 
patterns are easily produced in response to the details of 
how the samples are cut off at the ends. 

To test the response of the ground state to enhancing 
each of these elementary resonant loops, we slightly in- 
creased the exchange couplings along the bonds of such a 
loop at the center of a YC8 cylinder, and saw how much 
this enhanced the spin-spin correlations along the loop 
and elsewhere. It is the 8-site diamond loop that elicits 
the strongest response, as is shown in Fig. 5. The 6- 
site "perfect hexagon" loop that has received attention in 
previous papers^— shows a much smaller response, sug- 
gesting that this resonant loop is actually underweighted 
in the RVB ground state. 

One can not tile a kagome lattice with just these fa- 
vored resonant diamonds. However, a particular "dia- 
mond pattern" valence bond crystal, shown in Fig. 5, 
appears to be closely related to the spin liquid, and it is 
useful to think of the spin liquid as a melted state of this 
crystal. One must keep in mind that the actual over- 
lap of this state with the spin liquid is still small and 
the fluctuations are large. We have measured the corre- 
lation length for VBC correlations in this diamond pat- 
tern along cylinder YC8 and find that it is less than 1.5 
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FIG. 6: Ground state energy patterns for a YC10 (top) and 
a YC9-2 cylinder. The colors of the triangles, and their in- 
tensities, indicate the deviation of the sum of the spin-spin 
correlations on the three bonds forming the triangle from 3eo, 
where we take eo = —0.219. 



lattice spacings. Unlike the HVBC, this diamond VBC 
evolves smoothly into the spin liquid without any phase 
transition as one changes the strengths of the exchanges 
appropriately to favor it (we call this "pinning it"). For 
the even cylinders on which this diamond VBC does fit, 
this is very useful, as it allows a careful production of the 
spin liquid, by approaching it in a controlled and smooth 
fashion from this diamond VBC. 

This diamond pattern was essential in extracting our 
most accurate ground state energy estimate for YC12. 
Applying the pinning pattern of the lower panel of Fig. 
5 reduces the entanglement of the resulting state, allow- 
ing us to obtain more accurate energies. The pinning 
applies an equal number of positive and negative terms, 
so the energy dependence on the pinning coefficient rj 
has no linear term near rj = for the uniform spin liq- 
uid. After verifying this on several cylinders and de- 
termining that a nearly pure quadratic behavior held 
for rj < 0.02, we extrapolated E(rj) using the simple 
formula E(0) = 4/3£(0.01) - 1/3^(0.02). This proce- 
dure reduced our energy uncertainty for YC12 by a fac- 
tor of about 5. (The standard approach of simply ex- 
trapolating to zero truncation error at zero pinning gave 
-0.4375(14).) 

The infinitely long cylinders may be viewed as one- 
dimensional systems with a unit cell containing multi- 
ple spins. For those "even" cylinders that are com- 
patible with the diamond VBC, this unit cell contains 
an even number of spins (e.g., for YC8 it contains 12 
spins). In these cases the ground state of the infinite 
cylinder appears to be nondegenerate and gapped. But 
for the "odd" cylinders that are not compatible with 
the diamond VBC, the unit cell contains an odd num- 
ber of spins, and the Lieb-Schultz-Mattis theorem says 
the ground state must be degenerate.— We find that the 
ground states on these odd cylinders weakly break trans- 
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FIG. 7: Spin triplet (solid symbols) and singlet (hollow sym- 
bols) gaps for various cylinders with circumferences c. 



lational invariance, spontaneously doubling the unit cell, 
and this produces a pair of degenerate ground states, still 
with a gap to higher excited singlet states. The symme- 
try breaking is in a "striped" pattern that is shown in 
Fig. 6. For YC6 and YC10 the stripes run around the 
circumference, while for the other odd spiral cylinders 
the stripes are spirals. 



VI. GAPS 

To explore the low-lying excited states on our cylinders 
we use the following DMRG procedure: first target only 
one state, and sweep enough to obtain a high-accuracy 
ground state. Then restrict the range of bonds which 
are updated in the DMRG sweeps to the central half of 
the sample, and now target the two lowest-energy states, 
again sweeping to high accuracy, but keeping the end re- 
gions of the samples locally in the ground state. This re- 
stricted sweeping prevents the low-lying excitations from 
being bound to the ends of the sample. This technique is 
particularly important for obtaining the singlet gap. For 
the triplet gap, we can also apply magnetic fields on the 
ends to prevent any edge excitations from appearing and 
hiding the bulk gap. For the triplet, we can also target 
both states together, one with total S z = and the other 
with S z = 1, or run the two states separately. These dif- 
ferent approaches allowed for fairly independent checks 
on the results; in addition, we also varied the lengths and 
how quickly the number of states kept was increased. The 
results for these gaps are presented in Table I. Getting 
gaps is more demanding than getting ground state ener- 
gies, so our gap estimates do not go to as wide cylinders 
as do our ground state estimates. 

Our gap results for the singlet and triplet gaps differ 
in an interesting way from what is known about the gaps 
from exact diagonalizationj 25 i 29 

We believe that we have found good evidence for a 
nonzero singlet gap of about 0.04 or 0.05 in the 2D sys- 
tem. This is quite different from the exact diagonaliza- 
tion results, where there are many lower-lying singlets, 
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and the lowest singlet gap is only about 0.01 on the 
36-site torus^ As can be seen in Table I, the singlet 
gap is 0.050, within the errors, for the even XC8 and 
YC8 cylinders, and it remains near this value, although 
with much larger uncertainty, for the wider even XC12- 
2 cylinder. The odd cylinders, on the other hand, have 
smaller singlet gaps that increase as the circumference is 
increased, as shown in Fig. 7. The odd cylinders come in 
at least two "families": YC6 and YC10 are not spirals, 
while YC5-2, YC7-2 and YC9-2 are a series of spirals 
with increasing circumference. In each of these (short) 
sequences of odd cylinders the singlet gap increases as the 
circumference increases, supporting our conclusion that 
the singlet gap remains nonzero in the 2D limit. 

The triplet (spin) gap on the standard 36-site torus is 
0.164 from exact diagonalization^ While XC8 and YC8 
have gaps which are quite close to this, our results on 
other cylinders suggest that the two dimensional triplet 
gap is much smaller, as can be seen from Fig. 7. The 
triplet excitations appear to be composed of two spinons, 
but we cannot resolve whether or not the two spinons 
bind, although in some cylinders any binding must be 
very weak. This composite nature of the excitation seems 
to make the finite size effects and variation between the 
cylinders more pronounced. We do not yet understand 
the details of these effects. As in the singlet gap, the 
spiral odd cylinders have the smallest gaps, and the even 
the largest. Note that the triplet gap remains above the 
singlet gap in the systems we have studied, and we believe 
it remains nonzero in the 2D limit. 



VII. INDICATIONS OF Z 2 TOPOLOGICAL 
ORDER 

A nearest-neighbor RVB wavefunction is a linear com- 
bination of nearest-neighbor singlet dimers covering the 
kagome latticed For a kagome lattice wrapped on a 
cylinder, like we are studying, such dimer covers are in 
two topologically distinct classes, and dimer resonances 
on finite loops remain in the same topological sectori 21 i 22 
We can force our samples to be in one or the other of these 
two topological sectors by choosing how many spins to 
leave at each end. An example of this is apparent for 
the YC10 cylinder in Fig. 6, where we force the pattern 
by the number of spins we leave "sticking out" at each 
end. For odd cylinders, the two sectors are related by 
translation along the length of the cylinder, so are degen- 
erate. For even cylinders, on the other hand, the finite 
circumference c allows the two sectors to have different 
ground state energies. The difference in energy per site 
should vanish as ~ exp(— c/£), where £ is a correlation 
length. We can measure this splitting for YC8, where it 
is 0.00069(3) per site, but for the wider even cylinders we 
do not yet have an estimate, because it is more difficult 
to get a reliable ground state in the higher-energy topo- 
logical sector. It is also worth noting that on XC8 and 
YC8 the singlet gap above the higher-energy topological 



ground state is substantially smaller than the singlet gap 
above the overall ground state, which is presumably a 
factor in making the higher-energy sector more difficult 
to work with in DMRG. 

A domain wall along the cylinder between regions in 
the two different topological ground states is a spinon. 
For the odd cylinders, the degeneracy of the two topolog- 
ical sectors means the spinons are unconfined. However, 
two spinons might bind with a finite binding energy: it 
appears that this might be happening on cylinder YC10 
and thus it might also occur in the 2D limit. For the 
even cylinders, on the other hand, the spinons are con- 
fined by an effective potential that grows linearly with 
the distance, because the domain between them is in the 
higher-energy topological state. As a result, the excita- 
tion across the spin gap for the even cylinders must be a 
bound spinon pair. It remains to be determined if they 
remain bound in the 2D limit. 

We suspect that the lowest-lying "bulk" singlet exci- 
tations that we obtain are bound pairs of visons. It will 
be interesting to try to test this idea by attempting to 
pull apart such a vison pair or to try make a single vison 
near the end of a cylinder, but this remains for future 
research. 



VIII. CONCLUSIONS 

Using very accurate DMRG methods we have found 
a spin liquid ground state for the HKLM on cylinders 
with circumference up to 12 lattice spacings. The en- 
ergies are much lower than for the competing valence 
bond crystal state. The combination of low energies and 
small finite size effects due to short correlation lengths 
and nonzero singlet gaps, a new rigorous energy bound, 
and a simple picture for the nature of the spin liquid pro- 
vide compelling evidence that the infinite 2D system is 
a gapped RVB-like spin liquid. Much remains to be un- 
derstood concerning this phase, particularly the detailed 
structure, exchange statistics and dispersion relations of 
the various excitations. It should also be instructive to 
explore the phase diagram in the vicinity of this simple 
nearest-neighbor-only Heisenberg model by changing the 
Hamiltonian in various ways in order to find what other 
phases are nearby and perhaps to move "deeper" into 
this spin liquid where it might be easier to study. 
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